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ABSTRACT 

In this paper new real-time approaches for three-axis magnetometer sensor calibration are de- 
rived. These approaches rely on a conversion of the magnetometer-body and geomagnetic-reference 
vectors into an attitude independent observation by using scalar checking. The goal of the full cal- 
ibration problem involves the determination of the magnetometer bias vector, scale factors and 
non-orthogonality corrections. Although the actual solution to this full calibration problem in- 
volves the minimization of a quartie loss function, the problem can be converted into a quadratic 
loss function by a centering approximation. This leads to a simple batch linear least squares solu- 
tion. In this paper we develop alternative real-time algorithms based on both the extended Kalman 
filter and Unscented filter. With these real-time algorithms, a full magnetometer calibration can 
now be performed on-orbit during typical spacecraft mission-mode operations. Simulation results 
indicate that both algorithms provide accurate integer resolution in real time, but the Unscented fil- 
ter is more robust to large initial condition errors than the extended Kalman filter. The algorithms 
are also tested using actual data from the Transition Region and Coronal Explorer (TRACE). 

INTRODUCTION 

Three-axis magnetometers (TAMs) are widely used for onboard spacecraft operations. The basic 
concept behind these devices is a fairly simple one, involving a simple magnetic sensor coupled with 
an electronics unit to provide data in a digital format. These sensors are useful since they provide 
both the direction and magnitude of the magnetic field; they are lightweight, reliable, and have low 
power requirements, with no moving parts. 1 For these reasons most low-Earth orbiting spacecraft 
(below 1000 km) have TAMs as part of their basic sensor package. It is well known that a TAM 
can be used to determine a three-axis attitude when coupled with gyros or a dynamic model in 
an extended Kalman filter. 2 Attitude-knowledge accuracies of 1-2 degrees are common using this 
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approach, which can be improved by using well-calibrated sensors to achieve accuracies of 0. 1-0.5 
degrees. 3 An exciting new area of research involves using a TAM together with rate and sun sensor 
measurements for orbit (position) estimation. Accuracies on the order of 10 to 100 km can be 
achieved, which is within the position-knowledge requirements of many spacecraft. 4,5 These studies 
clearly show that an integrated magnetometer-based attitude/orbit estimation system can provide 
the necessary knowledge requirements of a spacecraft in a single package. 

A paramount issue to the attitude accuracy obtained using magnetometer measurements is the 
precision of the onboard calibration. The accuracy obtained using a TAM depends on a number 
of factors, including: biases, scale factors and non-orthogonality corrections. Scale factors and 
non-orthogonality corrections occur because the individual magnetometer axes are not orthonor- 
mal, typically due to thermal gradients within the magnetometer or to mechanical stress from the 
spacecraft. 6 We should note that even with perfectly calibrated (and noise-free!) magnetometers, 
inaccuracies still exist due to modelling uncertainty in the Earth’s geomagnetic field. These errors 
unfortunately are typically nonlinear and non-Gaussian in nature, which depend on many factors 
(see Ref. [7] for more details), and they ultimately limit the obtainable estimation performance us- 
ing a TAM. Magnetometer calibration is often accomplished using batch methods, where an entire 
set of data must be stored to determine the unknown parameters. This process is often repeated 
many times during the lifetime of a spacecraft in order to ensure the best possible precision obtained 
from magnetometer measurements. 

If an attitude is known accurately, then the magnetometer calibration problem is easy to solve. 
However, this is generally not the case. Fortunately, an attitude-independent scalar observation 
can be obtained using the norms of the body-measurement and geomagnetic-reference vectors. For 
the noise-free case, these norms are identical because the attitude matrix preserves the length of 
a vector. This process is also known as “scalar checking”. Unfortunately, even for the simpler 
magnetometer-bias determination problem, the loss function to be minimized is quartic in nature. 
The most common technique to overcome this difficulty has been proposed by Gambhir, who applies 
a “centering” approximation to yield a quadratic loss function that can be minimized using simple 
linear least squares. 8 Alonso and Shuster expand upon Gambhir’s approach by using a second step 
that employs the centered estimate as an initial value to an iterative Gauss-Newton method. Their 
algorithm, called “TWOSTEP”, 9 has been shown to perform well when other algorithms fail due 
to divergence problems. Furthermore, Alonso and Shuster have extended this approach to perform 
a complete calibration involving biases as well as scale factors and non-orthogonality corrections. 6 

One of the current goals of modern-day spacecraft is the ability to perform onboard and au- 
tonomous calibrations in real time without ground support. The disadvantage of the TWOSTEP 
algorithm is the requirement of a batch of data, which cannot be performed in real time. The 
centering approximation is linear and can be made into a sequential (real-time) algorithm, which is 
shown in this paper. However, the centering approach has been shown to be suboptimal for many 
realistic cases, 6 which is also confirmed by the results shown in this paper. The main objective 
of this paper is to present and compare two sequential algorithms that are suitable for real-time 
applications. The first algorithm uses an extended Kalman filter approach that is developed with 
commonly employed estimation techniques. The second algorithm uses an Unscented filter ap- 
proach that offers very good results for robust calibration when the initial conditions are poorly 
known. Experimental results using data obtained from the TRACE spacecraft show the validity of 
the new real-time algorithms to perform onboard and autonomous calibrations. 
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The organization of this paper is as follows: First, the attitude-independent observation model 
is summarized. Then, the centered linear algorithm for the solution of the complete magnetometer 
calibration problem is shown. This is extended to a sequential process, which can be used for 
real-time estimates. Next, a real-time algorithm based on the extended Kalman filter is derived, 
followed by an algorithm based on the Unscented filter. Finally, simulation and experimental results 
are shown. 

MEASUREMENT MODEL AND COVARIANCE 


In this section the TAM measurement model is summarized. Also, the attitude-independent 
observation as well as the covariance of the determined parameters are shown. More details on 
these concepts can be found in Ref. [6]. The magnetometer measurements can be modelled as 

B k = (I 3x3 + D)-\0 T A k H k + b + e k ), k = l,2,...,N (1) 


where B*. is the measurement of the magnetic field by the magnetometer at time H/~ is the 
corresponding value of the geomagnetic field with respect to an Earth-fixed coordinate system, A & 
is the unknown attitude matrix of the magnetometer with respect to the Earth-fixed coordinates, 
D is an unknown fully-populated symmetric matrix of scale factors (the diagonal elements) and 
non-orthogonality corrections (the off-diagonal elements), O is an orthogonal matrix (see Ref. [6] 
for a discussion on the physical connotations of this matrix), b is the bias vector, and e* is the 
measurement noise vector that is assumed to be a zero-mean Gaussian process with covariance E*. 
Also, I nxn is an n x n identity matrix. The goal of the full calibration problem is to estimate D 
and b. We first define the following quantities: 


E= 2D + D 2 

c = (-^3x3 + D)b 

Sk = [B\ k Bi Bl 2B lk B 2k 2B lk B 3k 2 B 2k B 3k ] 
E = [En B 22 E 33 E\ 2 Ei 3 £23 ] T 
An attitude-independent observation can be computed from 

z k = ||Bfc|| 2 — ||Hfc|| 2 

= L k e'-\\b(e ')\\ 2 + v k 

with 

L k = [2Bl -S k ] 


v k = 2[(I 3x3 + D)B k -b] T € k -\\€k\\ 2 


(2a) 

(2b) 

(2c) 

(2d) 


(3a) 

(3b) 


(4a) 

(4b) 

(4c) 


The effective measurement noise, v k , is approximately Gaussian with mean denoted by ^ k and 
variance denoted by a\, each given by 

Vk = E{v k } = - Tr(E fc ) (5a) 

a 2 k = E {vl} -& = 4[(/ 3 x3 + D) B fc - b] T E k [(I 3x3 + D) B k - b] + 2 (TV S 2 fc ) (5b) 
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Note that the measurement variance is a function of the unknown parameters. Equation (3a) is 
used to compute the actual effective measurement, while Eq. ( 3 b) represents the measurement 
model as a function of the unknown parameters E and c. 

Once E and c are determined, then a conversion from these parameters to D and b is required. 
This is accomplished by first decomposing E as 

E = UVU t (6) 

where U is an orthogonal matrix and V is a diagonal matrix with elements given by V\\, V22 and 
V33. We define W to be a diagonal matrix with elements given by 

Wjj = -1 + v'TTyJj, j = 1, 2 , 3 ( 7 ) 

Then, the matrix D is given by 

D = UWU t (8) 

Also, b can be computed using Eq. ( 2 b): 

b = (73x3 + -D) -1 c (9) 

Another quantity of interest is the transformation of the covariance of 6', denoted by to the 
covariance of 0, denoted by Poo , where 


9 = 


b 

D 


D = [-Dii -D22 D33 D\2 D 13 D23] ‘ 


(10a) 

(10b) 


The covariance Pg>o' is computed from a particular estimation algorithm. The covariance transfor- 
mation is given by 


where 


with 


Pee = 


d(b, D) 


[d(c, E) 


Pe'0' 


d(b, D) 


lT 


[d(c, E) 


\d(b, D)‘ 


r< 9 (c, E) 1 

-1 

(-^ 3 x 3 + D) 

M cD ( b) 


-1 

L«(c, E). 


kb, D)J 


06 x 3 

27 6 x 6 + M E d( D)_ 



2Du 0 

0 

2 D \2 

2 D 13 

0 

1 


0 

2D22 

0 

2 D \2 

0 

2D23 


0 

0 

2 T >33 

0 

2 D n 

2D23 

■d{P>) = 

D12 D\2 

0 

Dn + D22 

D 23 

D 13 


D\3 0 

D 13 

D 23 

Dn + D33 

D 13 


0 

D 23 

D 23 

D 13 

D\2 D22 + -D33 


( 11 ) 


( 12 ) 


( 13 ) 


4 


and 

M cD ( b) = 

Also, 0 mX n is an m x n matrix of zeros. 

CENTERED ALGORITHM 


b\ 0 0 62 63 0 

0 62 0 bi 0 63 

0 0 63 0 b\ 62 


(14) 


The measurement model in Eq. (3b) is clearly nonlinear in the unknown parameter vector O' . 
Therefore, linear least squares cannot be applied directly. However, it is possible to determine 
an approximate linear solution by applying a centering approach. First, the following weighted 
averages are defined: 



where 


1 =Y— 

* rr? 


k= 


Next, the following centered variables are defined: 


Lk = Lk — L 

Zk = Z k ~Z 

jl k = ilk - Ji 


(15a) 

(15b) 

(15c) 


(16) 


(17a) 

(17b) 

(17c) 


Then, an approximation for the optimal estimate of denoted by 0'*, is given by 6 



(18a) 

(18b) 


where Po^f is the covariance of O'*. Equations (6)-(9) can now be used to determine D and b. 
Also, Eq. (11) can be used to determine an approximation for the covariance of 0*, replacing P## 
with P$t 0 t. 


The approximate solution using Eq. (18) uses an entire batch of data. Since this solution is 
linear, then a sequential formulation can be derived that provides real-time estimates. A formal 
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derivation of this process can be found in Ref. [10]; we only present the final algorithm here. The 
sequential formulas for the averaged quantities in Eq. (15) are given by 

Lk+i = ~ 2 — —Z 2 (°fc+ l-kfc + H+i L k+i) (19a) 

°fc+ 1 + a k 

Zk+1 = _2 ~ -2 ipl+i^k + Vk+lZk+l) (19b) 

a k + 1 + a k 

Mfe+l = ~ 2 — —=2 + ^fe+lMfe+l) (19c) 

C7 k + 1 + 


where 

1 _ i 1 

-2 — -9 "I" 9 

4+1 * k 4+1 

Finally, the sequential formulas for & * and are given by 

= K k 0 k + ~~2 ( 4+1 ” 4 + 1 ) 

4+1 

Pe %+1 — K k P 9'6' k 

where the gain K k is defined by 


( 20 ) 


(21a) 

(21b) 


Kk = ^9x9 — Pe*e' k Lk+ 1 (4+i^^^fc+i + 4t+i) 4+i 


( 22 ) 


Note that only an inverse of a scalar quantity is required in the sequential process. Also, an 
approach for determining cr^ +1 involves using the previous estimate in Eq. (5b). The solutions 
for 0'x and P$*e f N can be shown to be equivalent to the batch solutions given by Eq. (18). The 
sequential process can be initialized using a small batch of data. A discussion on the observability 
of the centered estimate can be found in Ref. [11]. The TWOSTEP algorithm 9 uses the centered 
estimate as a starting value for a Gauss-Newton iteration (see Ref. [6] for more details), but this 
solution cannot be executed in real time. 


KALMAN FILTER FORMULATION 


In this section an extended Kalman filter (EKF) is derived to determine 0' in real time. A 
summary of the EKF equations can be found in Refs. [12] and [13]. Since the vector 0 f is constant, 
then the state model is given by 

i(t) - 0 (23) 

where x = 0'*, which is used to denote the optimal state estimate. The measurement model is 
given by 

z k = h(x fc ) + v k (24) 

where /i(xfc) = L k 0' k — ||bfc(0j[.)|| 2 from Eq. (3b). Since no process noise appears in the state 
model of Eq. (23), then the updated quantities (state and covariance) are given by their respective 
propagated quantities. The EKF equations then reduce down to 

4+1 = x fc -F K k [z k - h(x fc )] (25a) 

P k +i = [hx9-K k H k (x k )\P k (25b) 

K k = P k H k (xfc) [H k (yi k )P k Hl (x fc ) + <r|(x fc )] -1 (25c) 
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where P = Pe'e'- The state dependence of the measurement variance, cr|, is shown through Eq. (5b). 
The 1x9 matrix H (x) is the partial derivative of h(x) with respect to x. This quantity is given by 


H(x) = L- 


dllMflQIl 2 

ae 1 


(26) 


The partial derivatives of ||b(0')|| 2 are given by 6 


dllWH 2 

dE mn 


dCm 


= 2[(/ 3x3 +^r 1 c] m 


= -(2 - 6 mn ) [(/ 3 X3 + E)’ 1 c\ m [(/ 3x3 + E)- 1 c] n 


(27a) 

(27b) 


where [(/3x3 + £) 1 c] m denotes the m th element of (J 3x3 4- E) x c, and S mn is the Kronecker 
symbol. 


UNSCENTED FILTER FORMULATION 


In this section a new approach, developed by Julier, Uhlmann and Durrant- Whyte, 14 is shown 
as an alternative to the extended Kalman filter. This approach, which they called the Unscented 
filter (UF), typically involves more computations than the extended Kalman filter, but has several 
advantages, including: 1) the expected error is lower than the extended Kalman filter, 2) the 
new filter can be applied to non-differentiable functions, 3) the new filter avoids the derivation 
of Jacobian matrices, and 4) the new filter is valid to higher-order expansions than the standard 
extended Kalman filter. The Unscented filter works on the premise that with a fixed number of 
parameters it should be easier to approximate a Gaussian distribution than to approximate an 
arbitrary nonlinear function. The filter presented in Ref. [14] is derived for discrete-time nonlinear 
equations, where the system model is given by 

Xfc +1 = f(x fc , Wfc, k) (28a) 

z k = h(x fc , v fc , k) (28b) 

Note that a continuous-time model can always be written using Eq. (28a) through an appropriate 
numerical integration scheme. The vectors and are zero-mean Gaussian noise processes 
with covariances given by Q k and R k , respectively. The Kalman filter update equations axe first 
rewritten as 15 


Xfc = x fc + K k v k 

pt=p k -K k prK k 

where v k is the innovations process, given by 

v k = z k - z“ 

= z fc - h(x£ , k) 

The covariance of v k is defined by P^ v . The gain K k is computed by 

K k =pt(prr 1 ( 31 ) 

where P% z is the cross-correlation matrix between x^ and . 


(29a) 

(29b) 

(30) 
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The Unscented filter uses a different propagation than the standard extended Kalman filter. 
Given annxn covariance matrix P, a set of order n points can be generated from the columns 
(or rows) of the matrices ± VnP- The set of points is zero- mean, but if the distribution has mean 
/x, then simply adding p to each of the points yields a symmetric set of 2 n points having the 
desired mean and covariance . 14 Due to the symmetric nature of this set, its odd central moments 
are zero, so its first three moments are the same as the original Gaussian distribution. This is 
the foundation for the Unscented filter. A complete derivation of this filter is beyond the scope of 
the present paper, so only the final results are presented here. Various methods can be used to 
handle the process noise and measurement noise in the Unscented filter. One approach involves 
augmenting the covariance matrix with 



1 

p^ w 

p r 

II 

(Pr) r 

Qk 

p r 


X p k v ) T 

( p n T 

Rk _ 


(32) 


where P£ w is the correlation between the state error and process noise, Pg v is the correlation 
between the state error and measurement noise, and P™ v is the correlation between the process 
noise and measurement noise, which are all zero for most systems. Augmenting the covariance 
requires the computation of 2 (q + l) additional sigma points (where q is the dimension of w*. and 
l is the dimension of v*., which does not necessarily have to be the same dimension, m, as the 
output in this case), but the effects of the process and measurement noise in terms of the impact 
on the mean and covariance are introduced with the same order of accuracy as the uncertainty in 
the state. 


The general formulation for the propagation equations are given as follows. First, the following 
set of sigma points are computed: 

(Tk <— 2 L columns from ±7 y/Pj* (33a) 

xi( 0) = n (33b) 

X*(*) = o’* (*) + ** (33c) 


where L = n + q + 1 and xjj* is an augmented state defined by 






II 

X 

W* 

II 

a -a 
<X 

Ogxl 


_ v *. 


r-4 

X 

g 

o 


The parameter 7 is given by 

7 = Vl + x 

where the composite scaling parameter, A, is given by 

A = a 2 (L + K) - L 


(34) 


(35) 

(36) 


The constant a determines the spread of the sigma points and is usually set to a small positive 
value (e.g. 1 x 10 “ 4 < a < l ). 16 Also, the parameter k is usually given by k, — 3 — L. Efficient 
methods to compute the matrix square root can be found by using the Cholesky decomposition . 17 
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If an orthogonal matrix square root is used, then the sigma points lie along the eigenvectors of the 
covariance matrix. Note that there are a total of 2 L values for tr* (the positive and negative square 
roots). The transformed set of sigma points are evaluated for each of the points by 

Xk+i(i) = f (xjfe(i). Xk (0. k ) (37) 


where x£(i) is a vector of the first n elements of x£(i), and (*) is a vector of the next q elements 
of with 


X a k (i) = 


x ?(0 

Lx*(0j 


(38) 


where x'k(i) is a vector of the last l elements of which will be used to compute the output 

covariance. The following weights are now defined: 


WT- = J^r x (39a) 

W ™ = L TA + (1_a2+/?) (39b) 

^mea„ = ^cov = _l_ ) j = ^ 2 2 L (39c) 

where 0 is used to incorporate prior knowledge of the distribution (for Gaussian distributions 0 = 2 
is optimal). 


The predicted mean for the state estimate is calculated using a weighted sum of the points 
X%+ iW? which is given by 

2 L 

*w,=£ WT^xUiii) (40) 

1=0 

The predicted covariance is given by 


2 L 


n + 1 = E WE«M - %l WWO - %,f 


»=0 


The mean observation is given by 


where 


2 L 

t=0 


C*+l(») = Mxli-iW, Xjfe+i(0> fc+ 1) 

The output covariance is given by 

2 L 

PiU = E "T” Kwi (0 - %h! Kw-i(i) - 2 j +1 ] 3 

i= 0 


(41) 

(42) 

(43) 

(44) 


Then the innovations covariance is simply given by 


P >vv r>zz 

fc+l “ ^ 


Hi 


(45) 
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Finally the cross correlation matrix is determined using 

2 L 

Plh = E M '” v W - *i+J I 6 + 1 W - ^ +1 ] T (46) 

i = 0 

The filter gain is then computed using Eq. (31), and the state vector can now be updated using 
Eq. (29). Even though propagations on the order of 2n are required for the UF, the computations 
may be comparable to the EKF (especially if the continuous-time covariance equation needs to 
be integrated and a numerical Jacobian matrix is evaluated). Also, if the measurement noise, v*, 
appears linearly in the output (with l = m), then the augmented state can be reduced because the 
system state does not need to augmented with the measurement noise. In this case the covariance 
of the measurement error is simply added to the innovations covariance, with P%+i = P*+j + Pfc+ l* 
This can greatly reduce the computational requirements in the Unscented filter. Furthermore, a 
square root version of the UF is presented in Ref. [16] that avoids the need to re-factorize at each 
step. 

The implementation of the UF for real-time magnetometer calibration proceeds as follows. Two 
approaches in the UF are presented here. The first uses the same concept for the measurement-error 
variance as with the EKF, where the UF measurement-error variance is given by Rk = cr^(xfc). In 
this approach the augmented state vector is simply given by the state, so that a decomposition 
of a 9 x 9 matrix is only required, with Pg = P£. Since the state model estimate is given by 
x(t) = 0, then xj +1 = x£ and P k +j = P^, and Eqs. (40) and (41) are not needed (note: Qk = 0). 
This significantly reduces the computational requirements in the UF. Equation (43) is computed 
using the measurement model in Eq. (3b), with h(x*;, v*, k ) = Lk0' k — ||bjt(0j[.)|| 2 + v *. Also, the 
innovations covariance is given by P%+i = P k + 1 +Rk+ i* Hence, the only essential difference between 
the EKF and UF formulations is in the computation of the innovations covariance, where the EKF 
uses a first-order expansion to compute this quantity, while the UF uses a nonlinear transformation 
to compute this quantity. The second approach for the UF uses the measurement noise model of 
Eq. (3b) with an augmented vector given by the state and e in Eq. (4c). The augmented covariance 
is now given by 

'Pit 0 9x3 ‘ 

PS = (47) 

.0 3 x9 _ 

Therefore, a decomposition of a 12 x 12 matrix is now required. Again, since the state model 
estimate is given by x(t) = 0, then xjjT +1 = and P^ = P^", and Eqs. (40) and (41) are not 
needed. Also, the innovations covariance is computed using Eqs. (44) and (45) in this approach. 
In the strictest sense this approach is more optimal than the first approach because the effect of 
the nonlinear-appearing measurement noise is directly used in the UF. But, the computational 
requirements are vastly increased due to the decomposition of a higher dimensional augmented 
matrix. 

SIMULATION AND EXPERIMENTAL RESULTS 

In this section results of the EKF and UF formulations are shown using both simulated and 
experimental data. The simulated spacecraft is modelled after the Tropical Rainfall Measurement 
Mission (TRMM) spacecraft. This is an Earth-pointing spacecraft (spinning about its y-axis) 
in low-Earth orbit (currently near-circular at 402 km), with an inclination of 35°. The attitude 
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Table 1 Simulation Results 



Truth 

TWOSTEP 

Centered 

EKF 

UF 

b 1 

50 mG 

49.9620 ± 0.4320 

49.9700 ± 0.5525 

49.4029 ± 2.2625 

49.7289 ± 0.4328 

b 2 

30 mG 

29.8174 ± 0.4258 

29.5937 ± 5.5779 

19.5297 ±2.2044 

29.4657 ±0.3598 

b 3 

60 mG 

60.0331 ± 0.3165 

60.0456 ± 0.3165 

50.3248 ± 1.2340 

59.4561 ± 0.2331 

D n 

0.05 

0.0500 ± 0.0001 

0.0500 ± 0.0001 

0.0492 ± 0.0008 

0.0499 ± 0.0001 

D 22 

0.10 

0.0993 ± 0.0014 

0.0988 ± 0.0123 

0.0736 ± 0.0075 

0.0949 ± 0.0013 

D 33 

0.05 

0.0500 ±0.0001 

0.0500 ± 0.0001 

0.0481 ±0.0004 

0.0499 ± 0.0001 

D 12 

0.05 

0.0499 ± 0.0010 

0.0499 ±0.0011 

0.0313 ± 0.0051 

0.0486 ± 0.0010 

D 13 

0.05 

0.0499 ± 0.0001 

0.0499 ± 0.0001 

0.0440 ± 0.0006 

0.0495 ± 0.0002 

-D23 

0.05 

0.0501 ± 0.0007 

0.0501 ± 0.0007 

0.0287 ±0.0030 

0.0487 ±0.0006 


determination system consists of an Earth Sensor Assembly (ESA), a TAM, Digital Sun Sensors 
(DSS), Coarse Sun Sensors, and gyroscopic rate sensors. Currently, the ESA data is unreliable, 
so only the TAM and DSS are used as attitude sensors combined with gyros in an EKF (see 
Ref. [18] for more details). The geomagnetic field is simulated using a 10 th -order International 
Geomagnetic Reference Field model. 7 The magnetometer-body and geomagnetic-reference vectors 
for the simulated runs each have a magnitude of about 500 milliGauss (mG). The measurement 
noise is assumed to be white and Gaussian, and the covariance is taken to be isotropic with a 
standard deviation of 0.5 mG. The measurements are sampled every 10 seconds over an 8-hour 
span. The true values for the bias b and elements of the D matrix are shown in Table 1. Large 
values for the biases are used to test the robustness of the EKF and UF algorithms. 

Thirty runs have been executed, which provide a Monte-Carlo type simulation. Shown in Table 
1 Eire the averaged batch solutions given by the TWOSTEP and centered algorithms, each with 
their maximum deviations obtained. Since the TWOSTEP approach is the most rigorous, all 
comparisons are made with respect to this algorithm. The centered algorithm does a fairly good 
job at estimating all parameters, with the exception of 62- This parameter corresponds to the least 
observable variable, which results in a wide variation from the averaged value. 


The EKF and UF are both executed at time t = 0 using initial conditions of zeros for all states. 
The initial covariance matrix is diagonal, given by 


Po = 


500/3x3 

06x3 


03x6 

0.001/6x6. 


(48) 


This assumes a 3cr bound on the initial bias estimates to be about 70 mG and a 3 o bound on the 
initial estimates for the elements of the D matrix to be about 0.1. The parameters used in the UF 
are a = 0.1, {3 — 2, k = 3 — L, and L — 9 when the 9x9 matrix decomposition approach is used 
(L = 12 for the augmented approach). The two UF approaches (one uses the current state estimate 
in the measurement-error variance calculation, which requires a 9 x 9 matrix decomposition, while 
the other appends the state vector to include the measurement noise, which requires a 12 x 12 matrix 
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Figure 1 EKF and UF Estimates and 3 a Bounds for the Parameter 63 

decomposition) give identical results. Therefore, no apparent advantages to using the augmented 
approach in the UF is seen. The EKF and UF solutions at the final time are shown in Table 1. 
The EKF does not converge to the correct solution for many of the parameters, while the UF gave 
results that are just as good as the TWOSTEP solutions. Also, the maximum deviations for the 
UF are much smaller than for the EKF and centered algorithm. 

Figure l(i) shows the EKF and UF estimates for the parameter 63 . The EKF does not converge 
to the correct solution during the 8 -hour simulated run. This is due to the fact that the first-order 
approximation in the EKF does not adequately capture the large initial errors. The biggest concern 
with the EKF solutions is the confidence of the results dictated by the 3 a bounds, with 63 shown 
in Figure l(ii). In fact, this plot shows that the EKF is performing better than the UF. This 
can certainly provide some misleading results. However, unlike the EKF, the maximum deviations 
associated with the UF shown in Table 1 are within the 3 cr bounds for all the parameters. This 
indicates that the UF is performing in an optimal manner. But, the UF algorithm comes with 
a computational cost, m ainl y due to the covariance decomposition. Our experience has shown 
that the UF algorithm is about 2 times slower than the EKF algorithm. Still, the performance 
enhancements of the UF over the EKF may outweigh the increased computational costs. 

Next, experimental results from the TRACE spacecraft are shown. This is an Sun-synchronous 
spacecraft in low-Earth orbit (currently near-circular at 402 km). The attitude determination 
system consists of a TAM, DSS and gyroscopic rate sensors. The data collected for the spacecraft 
is given during an inertial pointing mode. As discussed previously the errors associated with the 
geomagnetic field model are typically nonlinear and non-Gaussian in nature. This violates the 
assumptions for all the estimators shown in this paper. We still assume that the measurement 
noise is white and Gaussian, but the standard deviation is now increased to a value of 3 mG, which 
bounds the errors in a practical sense. The measurements are sampled every 3 seconds over a 
6 -hour span. Figure 2 (i) shows the magnetic field reference for all axes, and Figure 2 (ii) shows the 
actual TAM measurements. 
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Figure 2 Reference and Measurements from TRACE, and Bias Estimates and 3cr Bounds 


The EKF and UF are both executed at time t = 0 using initial conditions of zeros for all states. 
The initial covariance matrix is diagonal, given by 


Po = 


fl0/3x3 


06x3 


03x6 


0.001/6x6 


(49) 


This assumes a 3cr bound on the initial bias estimates to be about 10 mG and a 3cr bound on the 
initial estimates for the elements of the D matrix to be about 0.1. The parameters used in the 
UF are a = 0.1, /? = 2, k — 3 — L, and L = 9 when the 9x9 matrix decomposition approach 
is used (L = 12 for the augmented approach). Once again no advantages to using the augmented 
approach in the UF are seen. For the experimental data the solutions obtained using TWOSTEP, 
the EKF and UF algorithms are nearly identical. This is most likely due to the well-behaved nature 
of the data (i.e., the calibration errors are small). However, the sequential centered algorithm gave 
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slightly different results, which indicates that this approach is not truly optimal, even when small 


errors axe present. The centered algorithm final results are given by 

b*= [1.4007 -8.7350 -3.7927 ] T (50a) 

D* = [0.0086 0.0437 0.0065 0.0006 0.0035 -0.0120] T (50b) 

The TWOSTEP, EKF and UF final results are given by 

b*= [1.6056 -8.4140 -4.6123 ] T (51a) 

D* = [0.0123 0.0181 0.0040 -0.0005 0.0038 -0.0019] T (51b) 


These results indicate that the EKF and UF algorithms should be used over the sequential centered 
algorithm for real-time calibration purposes. Figure 2(iii) shows EKF estimates for the bias vector 
b. Another advantage of a real-time approach is the convergence properties of the particular 
estimator. From Figure 2(iii) good convergence is seen for all the parameters, which indicates that 
the calibration parameters are well behaved (i.e., truly constant in a practical sense). Figure 2(iv) 
shows the 3cr bounds for the bias estimates. This at least qualitatively indicates that good parameter 
estimates are achieved since these bounds are fairly small compared to the TAM measurements. 
Similar results are obtained for the D matrix parameters. 

An investigation of the residuals between the norm of the estimated vector, using the calibrated 
parameters, and the geomagnetic-reference vector is useful to check the consistency of the results. A 
plot of this residual is shown in Figure 3. From this plot the residuals are slightly biased by about 2 
mG and have some non-Gaussian components. This clearly shows the nonlinear and non-Gaussian 
nature of the geomagnetic field model. A spectrum analysis shows the presence of sinusoidal motions 
with periods equivalent to the orbital period (~ 90 min) and higher-order harmonics (see Ref. [19] 
for a model of this process). Still, the residuals are small enough to indicate that good results are 
obtained from the real-time calibration algorithms developed in this paper. 
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CONCLUSIONS 


In this paper two new real-time algorithms based on the extended Kalman filter and Unscented 
filter were developed for the calibration of three-axis magnetometers. Also, a sequential algorithm 
was shown that was based on a linear least squares approach using a centering approximation. Two 
different approaches have also been shown in the Unscented filter design. One uses the current state 
estimate in the measurement-error variance calculation, which requires a 9x9 matrix decomposition, 
while the other appends the state vector to include the nonlinear-appearing measurement noise, 
which requires a 12 x 12 matrix decomposition. Simulation and experimental results indicated 
that for this problem, both Unscented filter approaches gave identical results. So appending the 
state vector is not required, which significantly reduces the computational requirements in the 
Unscented filter. Simulation test cases also showed that the performance of the Unscented filter 
is significantly better than the standard extended Kalman filter and sequential centered algorithm 
for large initialization errors. Experimental results indicated that the EKF and UF algorithms 
should be used over the sequential centered algorithm for real-time calibration purposes, even 
when small errors are present. However, the Unscented filter algorithm is recommended for actual 
implementation when computational requirements are not burdensome. 
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